Spurious diffusion in particle simulations of the Kolmogorov flow 
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Abstract 



Particle simulations of the Kolmogorov flow are analyzed by the Landau- 
Lifshitz fluctuating hydrodynamics. It is shown that a spurious diffusion 
of the center of mass corrupts the statistical properties of the flow. The 
analytical expression for the corresponding diffusion coefficient is derived. 
Some fifty years ago, Kolmogorov introduced a simple hydrodynamic model, that displays 
various instabilities, including a transition to turbulence |JJ. The richness of this model, 
known as the Kolmogorov flow, combined with its simplicity, has attracted a great number 
of both theoretical and numerical work [0-|5| . In particular, particle simulations, based on 
Lattice Boltzmann || or Lattice gas automata [jTJ, have been used to study the statistical 
properties of the flow in high Reynolds number regime iPJsHTofl . The purpose of this letter 
is to point out a subtle problem concerning the nonequilibrium fluctuations that appear in 
this model. We show that the center of mass of the system undergoes a spurious diffusion 
that corrupts the statistical properties of the flow. 

The Kolmogorov flow is an isothermal fluid confined in a rectangular box L x x L y , 
{0 < x < L x , < y < L y }, with periodic boundary conditions in both directions. The flow 
is maintained through an external force field of the form 



where l x is the unit vector in the x direction. The relevant parameters are the strength of 
the force field F , the wave number n of the forcing, and the aspect ratio a r = L x /L y . 

For small enough F , the flow follows basically the external field and the stationary 
velocity profile is readily found to be 
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where 77 is the shear viscosity coefficient. However, upon increasing Fq, this stationary 
state becomes unstable giving rise to rotating convective patterns 0,0]. Other instabilities 
of increasing complexity appear for larger values of Fq, culminating in a turbulent - like 
behavior ||||. 

To study the statistical properties of this system, we turn to a description in terms of 
the fluctuating hydrodynamic equations: 

§f = -V ■ (,v) (4) 
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where p is the mass density, P the hydrostatic pressure and cr the two-dimensional stress 
tensor: 

= -V (f^+f^-^ v - v ) (6) 

S is a random tensor whose elements {Sij} are Gaussian white noises with zero mean and 
covariances given by JTT 
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where ks and T stand for the Boltzmann constant and the (uniform) temperature, respec- 
tively. For simplicity, we shall assume that the shear and bulk viscosity coefficients, 77 and 
C, are state independent, i.e., they are constant. 

When imposing a force field, one has to keep in mind that both in microscopic simulations 
and in real systems, the fluid is made out of individual particles. Hence one cannot impose 
a bulk force, but rather an acceleration field acting on the particles. Since the density of 
particles is fluctuating, we conclude that the external field in the momentum equation (|5|) 
is also a fluctuating quantity: 

F ext = p(x, y) a sin (2 7r n y/ L y ) l x , (8) 



where ao is the amplitude of the imposed acceleration field. Furthermore, since the external 
field in the Kolmogorov flow is space-dependent, the force acting on a particle depends on 
its exact position so that the total force (in the x-direction) F(t) will also be fluctuating, 
even though the total number of particles is conserved. As a result, the center of mass linear 
momentum, denoted by J x (t), undergoes a stochastic motion driven by a scalar force F(t): 

dJjJ)_ = rL. ^ rL« rf ( j sm(2nny/L y ) = F(t). (9) 
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In strictly subsonic regimes the flow behaves essentially as an incompressible fluid so 
that the average density is uniform in space, < p >= po- It then follows from eq. @ 
that < F(t) >= 0. To find the force correlation function < F(t)F(t') >, we consider the 
spatial average of the hydrodynamical equations [5]) over the x direction and notice that 
corresponding spatially averaged density, p{y,t) = j x dx p(x,y,t), and y component of 
the velocity, v(y,t) = J Lx dxv y (x,y,t), are not affected by the external constraints, i.e. 
they assume their equilibrium form. In particular, in the stationary regime one has that 
< p >= Po and < v >= are independent of the value of ao- To study the fluctuations 
around this state, we introduce the deviations Sp(y,t) = p(y,t) — p , 5v(y,t) = v(y,t) and 
5P(y,t) = P{y,t) — < P >, that obey the following linearized equations : 

d8p dSv , . 
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with 
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To close these equations, we need to specify the equation of state. Since the fluid is isother- 
mal, we simply set 

5P = c 2 s Sp, (13) 
where c s is the isothermal sound speed. 



The stochastic differential evolution equation for the fluctuating force F(t) now follows 
easily by multiplication of (|T0| ) and (|TTJ) by sin (2 ir n y/L y ) and cos {2tx ny j L y ), respectively, 
followed by integration over y. One obtains: 

d 2 F(t) V + ( W dF(t) Wcg _ 
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where ip(t) is a white Gaussian noise, with zero mean and variance given by: 

k B T 2 /2vrn x 
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We conclude that is a Gaussian non-Markovian process. The exact form of the force 
correlation is easily obtained from ( [H] ) and flT5f), but the final expression is rather lengthy. 
On the other hand, the validity of hydrodynamics can only be guaranteed if the parameter 

v + C 
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remains small |fL2|| . Accordingly, to dominant order in e, the force correlation reads: 
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< F(t)F(0) >= p 2 Q al B ° exp {- 4ttV T s t/Lj} cos(27m c s t/L y ) ; t > 0, (17) 

where r s = (r/ + ()/2p represents the (two dimensional) sound damping coefficient, N is 
the total number of particles and m their individual mass. 

Turning to J x {t), which is nothing but the time integral of F(t), we conclude that it is a 
Gaussian stochastic process with zero average and second moment (again to dominant order 
in e): 

+ (27m/ L y ) 1 — exp (-4:7i 2 n 2 T s t/L 2 y ) cos(2irnc s t/L y )\ J (18) 

As announced, J x (t) diffuses in time (in the momentum space) with a long-time diffusion 
coefficient given by: 
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It is important to notice that in real macroscopic systems the very existence of the 
center of mass diffusion remains questionable for the following reason. Periodic boundary 
conditions are one of the basic simplifying features of the Kolmogorov flow. This is fine 
for the representation of a system of infinite extend, consisting of periodically repeated 
Kolmogorov units, as long as only macroscopic properties are at stake. However, when the 
fluctuations are under investigation, one has to realize that the periodic boundary conditions 
imply a perfect correlation of the fluctuating forces in the different units. This is obviously 
unphysical (except for the academic situation of a system defined on a torus). In any case, 
the diffusion coefficient D is negligibly small in macroscopic systems (typically of the order 
of 10~ 38 A;(? 2 /m 2 s 3 ) and is essentially unobservable. 

The situation is entirely different in microscopic simulations where the total number of 
particles N barely exceeds 10 5 . To estimate the importance of the center of mass diffusion, 
and the corresponding contamination of statistical properties of the system in numerical 
simulations, we first note that the ratio ks T /mc^ is of the order of unity. We next observe 
that there is a minimum run time for simulations, namely the hydrodynamic relaxation time 
Th pa L x L y /T s . Typical running times are several such r^. It then follows from eq. fllBD that 
for large t (i.e. t > r h ) the center of mass velocity fluctuation, < J x > /pi, is about: 

<vl{t)>= < Jl f > ^ A-L, (20) 

where n = N / L x L y is the number density. This quantity has to be compared with the 
spatial average of the mean square flow velocity u^, which is of the order of (see eq. 
(^)). The relative importance of the center of mass diffusion can thus be estimated by the 
square root of the ratio < v x {t) > /u^, hereafter denoted by //(£). Using the explicit form 
of uq, eq. (0), one finds that for t > r^, 

where a r = L x /L y is the aspect ratio. 

As a first example, we consider a two dimensional Boltzmann gas for which there exists 
an efficient algorithm, proposed two decades ago by Bird, that is about 3 orders of magnitude 



( < vl{t) > /4J 1/2 - (2/no) V V (t/r,) 1/2 (21) 



faster than the corresponding traditional molecular dynamic simulation |fL3]| . A typical case 
is a system involving 20 000 hard disks of diameter d, with L x x L y = 2000 x 1000 d 2 (i.e. 
a r = 2 and no = 10~ 2 particles per d 2 ), n = 2, c s ~ 1 and 77 ~ 0.3 (in system units, where 
lengths, masses and velocities are scaled by the disk diameter d, the particle mass m and 



the thermal velocity, y^To/m, respectively). It then follows from eq. (|ZID that after only 
one relaxation time, /i(r^) ~ 7 x 10 -2 which is certainly not negligible, all the more so since 
typical running times are 10 to 100 times larger than r^. 

One way to avoid this problem is to increase the number of particles, while keeping the 
number density uq = 10~ 2 particles per d 2 , since then the Bird algorithm is applicable. 
However, to reach reasonably small values of fi, like for instance /^(t/J ~ 10~ 4 , one has 
to consider simulations involving over 13 millions of particles. Such simulations require a 
prohibitively long running time with present day computers. 

The only other alternative is to increase the number density as well. For a given number 
of particles, the best strategy is to choose n so that the Reynolds number is as high as 



possible, since this is precisely one of the main objectives of numerical simulations |jT3 
In the case of subsonic hard disk flows, the appropriate number density turns out to be 
about n = 0.27 particles per d 2 ||15|. For a system containing half a million of particles, 



L x x L y = 960 x 1920 d 2 , c s « 1.6 and 77 ps 0.4. The function ^{j h ) is then about 4 x 10~ 4 , 
which is quite satisfactory. However, a number density of n Q « 0.27 corresponds to a 
moderately dense Enskog gas for which the Bird algorithm is no longer applicable |]16" 



Instead, one has to use the traditional hard disk molecular dynamics method which, as 
mentioned before, is about 3 orders of magnitude slower than the corresponding dilute gas 
simulation. Furthermore, the collision frequency grows linearly with the number density, 
which further increases the running time by at least another order of magnitude. Under 
these conditions, pursuing the simulation for a single relaxation time 77, is about the best 
one can achieve with present day computer performances. Although such a relatively short 
simulation might be satisfactory to study the average properties of the system, it is certainly 
not enough to extract the associated fluctuation spectrum. 



The above discussion highlights the usefulness of lattice-particle simulations for the study 
of the relatively high Reynolds number flows. But these model simulations have their own 
limitations. Because the motion of particles takes place within a restricted geometry (4 or 
6 linear directions), with the corresponding restricted number of velocities, reaching local 
equilibrium requires now many more collisions than in the case of hard disk dynamics |T7] . 
As far as macroscopic properties of the system are concerned, this is only a minor problem, 
since lattice-particle simulations typically run seven orders of magnitude faster than hard 
disk molecular dynamics. The major drawback however is that such a long time simulation 
inevitably increases the effect of the center of mass diffusion reported here. In fact, the 
spurious diffusion has been noticed very recently by Boon et al. |L8| in a study of the so- 
called "turbulent diffusion" in Kolmogorov flow. 

In conclusion, while spurious diffusion of the center of mass in the Kolmogorov flow 
does not affect the average macroscopic behavior of the system, it does corrupt the other 
statistical properties, and to a significant degree under conditions that are typical for many 
microscopic simulations. The best way to avoid this problem is to include in the simulation 
algorithm an ad-hoc mechanism that prevents the center of mass momentum fluctuations. 



This can be accomplished rather easily in lattice-particle simulations [I5|, but its counterpart 
in molecular dynamic simulations is less obvious. 



I. ACKNOWLEDGMENTS 

We thank G. Nicolis and J. W. Turner for stimulating discussions. We acknowledge 
support from the IUAP, Prime Minister's Office of the Belgian Government, and one of 
us (CVdB) also from the FWO Vlaanderen. This work has been partly supported by the 
European Commission DG 12 (Advanced Research Meetings and Studies). 



8 



REFERENCES 

[1] A.M. Obukhov, Russ. Math. Survey 38,113 (1983). 

[2] L. Meshalkin, and Ya.G. Sinai, J.Appl.Math.Mech. 25, 1140 (1961) ; J.S.A. Green, 
J.Fluid Mech. 62, 273 (1974). 

[3] A. A. Nepomniaschichii, J. Appl. Math. Mech. 40, 836 (1976). 

[4] E. N. Lorenz, J. Fluid Mech., 55, 545 (1972). 

[5] Z. S. She, Phys. Lett. A 124, 161 (1987) ; ibid. Ph.D. Thesis, Univ. Paris VII (1987) ; 
G. I. Sivashinsky, Physica D 17, 243 (1985). 

[6] R. Benzi, S. Succi and M. Vergassola, Phys. Reo. 222, 147 (1992). 

[7] U. Frisch, B. Hasslacher and Y. Pomeau, Phys. Rev. Lett. 56, 1505 (1986) ; G. Doolen, 
Ed., Lattice Gas Methods for PDE's, Physica D 47, 1 (1991). 

[8] M. Henon and H. Scholl, Phys. Rev. A 43, 5365 (1991). 

[9] R. Benzi, S. Succi, J. Stat. Phys. 56, 69 (1989) ; U. Frisch, B. Legras, B. Villone, 
Physica D 94, 36 (1996). 

[10] E. Vanden Eijnden, D. Hanon and J. P. Boon, in Advances in Turbulence, U. Frisch, 
ed., VII, 183 (1998). 

[11] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Pergamon Press, Oxford (1984). 

[12] See for instance W. E. Alley and B. J. Alder, Phys. Rev. A 27, 3158 (1983) ; W. E. 
Alley, B. J. Alder and S. Yip, Phys. Rev. A 27, 3174 (1983). 

[13] G. A. Bird, Molecular Gas Dynamics, Clarendon, Oxford (1976) ; ibid Gas Dynamics 
and the Direct Simulation of Gas Flows, Clarendon, Oxford (1994). 

[14] A Bird algorithm with adjustable transport coefficients has been proposed in : F. Baras, 
M. Malek Mansour. A. Garcia and M. Mareschal, J. Comp. Phys. 119, 94 (1994). 

9 



[15] A. Puhl, M. Malek Mansour and M. Mareschal, Phys. Rev. A 40, 1999 (1989). 

[16] An extension of the Bird algorithm dealing with dense systems has been reported in : 
F. Alexander, A. Garcia and B. Alder, Phys. Rev. Lett. 74 5212 (1995). 

[17] D. H. Rothman and S. Zaleski, Lattice-Gas Cellular Automata, Cambridge University 
Press (1997). 

[18] J. P. Boon, E. Vanden Eijnden and D. Hanon (to appear in Chaos, Solitons and Fractals). 
[19] D. Hanon (private communication). 



10 



